Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

6.3.1.2 Number of samples used

[1] 27
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.3 Richness


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.972
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.973
Podarcis_muralis            0.97302                 
Df SumOfSqs R2 F Pr(>F)
species 1 1.542719 0.2095041 6.625717 0.001
Residual 25 5.820951 0.7904959 NA NA
Total 26 7.363669 1.0000000 NA NA

6.3.1.4 Neutral


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.941
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.944
Podarcis_muralis            0.94763                 
Df SumOfSqs R2 F Pr(>F)
species 1 1.918266 0.2608511 8.822682 0.001
Residual 25 5.435610 0.7391489 NA NA
Total 26 7.353876 1.0000000 NA NA

6.3.1.5 Phylogenetic


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.132
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                               0.13
Podarcis_muralis            0.12705                 
Df SumOfSqs R2 F Pr(>F)
species 1 0.3218613 0.2162815 6.899207 0.001
Residual 25 1.1662981 0.7837185 NA NA
Total 26 1.4881594 1.0000000 NA NA

6.3.1.6 Functional


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.018367 0.018367 1.5597    999  0.205
Residuals 25 0.294402 0.011776                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.199
Podarcis_muralis            0.22328                 
Df SumOfSqs R2 F Pr(>F)
species 1 0.0858578 0.172879 5.225323 0.057
Residual 25 0.4107775 0.827121 NA NA
Total 26 0.4966352 1.000000 NA NA
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")

6.3.2.1 Number of samples used

[1] 27
beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.2 Richness


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.11796 0.117959 12.963    999  0.003 **
Residuals 25 0.22748 0.009099                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  0.0013711        
Df SumOfSqs R2 F Pr(>F)
Population 1 1.639807 0.179834 5.481634 0.001
Residual 25 7.478640 0.820166 NA NA
Total 26 9.118447 1.000000 NA NA

6.3.2.3 Neutral


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07844 0.078443 5.2384    999  0.032 *
Residuals 25 0.37437 0.014975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.029
Hot_dry  0.030815        
Df SumOfSqs R2 F Pr(>F)
Population 1 1.947003 0.2306127 7.493387 0.001
Residual 25 6.495736 0.7693873 NA NA
Total 26 8.442739 1.0000000 NA NA

6.3.2.4 Phylogenetic


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.06739 0.067395 2.9532    999  0.092 .
Residuals 25 0.57052 0.022821                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.088
Hot_dry  0.098068        
Df SumOfSqs R2 F Pr(>F)
Population 1 0.2441653 0.1224638 3.488854 0.019
Residual 25 1.7496100 0.8775362 NA NA
Total 26 1.9937754 1.0000000 NA NA

6.3.2.5 Functional


Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02496 0.024955 0.6729    999   0.44
Residuals 25 0.92714 0.037085                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.432
Hot_dry   0.41979        
Df SumOfSqs R2 F Pr(>F)
Population 1 0.0279454 0.0248037 0.6358634 0.444
Residual 25 1.0987171 0.9751963 NA NA
Total 26 1.1266624 1.0000000 NA NA
beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

6.3.3.1 Number of samples used

[1] 54
beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)
6.3.3.1.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05014 0.050145 6.2252    999  0.027 *
Residuals 52 0.41886 0.008055                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                          0.02
1_Acclimation 0.015808              
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.6172653 0.0360987 2.320618 0.001
species 1 2.8279677 0.1653842 10.631785 0.004
time_point:species 1 0.3545578 0.0207351 1.332965 0.054
Residual 50 13.2995905 0.7777820 NA NA
Total 53 17.0993812 1.0000000 NA NA
6.3.3.1.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0199 0.0199035 2.1213    999  0.149
Residuals 52 0.4879 0.0093827                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.156
1_Acclimation 0.15128              
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.9050519 0.0541893 3.792749 0.001
species 1 3.3236300 0.1989999 13.928143 0.001
time_point:species 1 0.5416391 0.0324302 2.269815 0.004
Residual 50 11.9313461 0.7143805 NA NA
Total 53 16.7016671 1.0000000 NA NA
6.3.3.1.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01334 0.013340 0.6524    999  0.407
Residuals 52 1.06332 0.020449                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.403
1_Acclimation 0.42294              
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.2890434 0.0766494 4.956318 0.002
species 1 0.3508889 0.0930498 6.016803 0.001
time_point:species 1 0.2151377 0.0570509 3.689034 0.008
Residual 50 2.9159082 0.7732498 NA NA
Total 53 3.7709782 1.0000000 NA NA
6.3.3.1.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0123 0.012300 0.4817    999  0.504
Residuals 52 1.3277 0.025533                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                         0.51
1_Acclimation 0.49073              
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.0448774 0.0269021 1.4865056 0.248
species 1 0.0973005 0.0583275 3.2229509 0.346
time_point:species 1 0.0165026 0.0098926 0.5466273 0.397
Residual 50 1.5094945 0.9048777 NA NA
Total 53 1.6681751 1.0000000 NA NA
beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Do the antibiotics work?

6.3.4.1 Antibiotics

treat1 <- meta  %>%
  filter(time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat1))]
identical(sort(colnames(treat1.counts)),sort(as.character(rownames(treat1))))

treat1_nmds <- sample_metadata %>%
  filter(time_point == "2_Antibiotics")

6.3.4.2 Number of samples used

[1] 23
beta_div_richness_treat1<-hillpair(data=treat1.counts, q=0)
beta_div_neutral_treat1<-hillpair(data=treat1.counts, q=1)
beta_div_phylo_treat1<-hillpair(data=treat1.counts, q=1, tree=genome_tree)
beta_div_func_treat1<-hillpair(data=treat1.counts, q=1, dist=dist)
6.3.4.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.015319 0.0153186 6.8764    999  0.015 *
Residuals 21 0.046782 0.0022277                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.009
Hot_dry  0.015919        
Df SumOfSqs R2 F Pr(>F)
Population 1 1.356644 0.1527052 3.784762 0.001
Residual 21 7.527429 0.8472948 NA NA
Total 22 8.884073 1.0000000 NA NA
6.3.4.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.030536 0.0305358 3.8593    999  0.057 .
Residuals 21 0.166158 0.0079123                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.067
Hot_dry  0.062842        
Df SumOfSqs R2 F Pr(>F)
Population 1 1.785669 0.2085055 5.532084 0.001
Residual 21 6.778468 0.7914945 NA NA
Total 22 8.564137 1.0000000 NA NA
6.3.4.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.012041 0.012041 0.9898    999  0.359
Residuals 21 0.255459 0.012165                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.345
Hot_dry   0.33111        
Df SumOfSqs R2 F Pr(>F)
Population 1 0.8963254 0.1888758 4.889993 0.001
Residual 21 3.8492558 0.8111242 NA NA
Total 22 4.7455811 1.0000000 NA NA
6.3.4.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01802 0.018021 0.4386    999  0.484
Residuals 21 0.86280 0.041086                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.485
Hot_dry   0.51499        
Df SumOfSqs R2 F Pr(>F)
Population 1 0.0184663 0.0098404 0.2087022 0.732
Residual 21 1.8581156 0.9901596 NA NA
Total 22 1.8765819 1.0000000 NA NA
beta_richness_nmds_treat1 <- beta_div_richness_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat1 <- beta_div_neutral_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat1 <- beta_div_phylo_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat1 <- beta_div_func_treat1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat1_nmds, by = join_by(sample == Tube_code))

6.3.4.3 Acclimation vs antibiotics

treat <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

6.3.4.4 Number of samples used

[1] 50
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)
6.3.4.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.025318 0.0253178 6.021    999  0.014 *
Residuals 48 0.201837 0.0042049                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.012
2_Antibiotics      0.017817              
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.8885838 0.0949462 5.789315 0.001
Population 1 2.1171094 0.1064350 6.489843 0.001
time_point:Population 1 0.8793415 0.0442078 2.695557 0.002
Residual 46 15.0060684 0.7544111 NA NA
Total 49 19.8911031 1.0000000 NA NA
6.3.4.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.039587 0.039587 6.8387    999  0.007 **
Residuals 48 0.277854 0.005789                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                        0.01
2_Antibiotics      0.011886              
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.0241808 0.1063620 7.014531 0.001
Population 1 2.8531033 0.1499183 9.887052 0.001
time_point:Population 1 0.8795688 0.0462175 3.048030 0.003
Residual 46 13.2742044 0.6975022 NA NA
Total 49 19.0310573 1.0000000 NA NA
6.3.4.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.58372 0.58372 35.413    999  0.001 ***
Residuals 48 0.79119 0.01648                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.001
2_Antibiotics    2.9795e-07              
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.8065206 0.2113909 14.842282 0.001
Population 1 0.7903334 0.0924813 6.493340 0.001
time_point:Population 1 0.3501572 0.0409738 2.876874 0.043
Residual 46 5.5988658 0.6551540 NA NA
Total 49 8.5458771 1.0000000 NA NA
6.3.4.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18591 0.185914 5.0679    999  0.033 *
Residuals 48 1.76088 0.036685                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.027
2_Antibiotics      0.028989              
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.8020952 0.3750193 28.0355332 0.001
Population 1 0.0031247 0.0006503 0.0486115 0.001
time_point:Population 1 0.0432870 0.0090081 0.6734238 0.482
Residual 46 2.9568327 0.6153223 NA NA
Total 49 4.8053396 1.0000000 NA NA
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Does the FMT work?

6.3.5.1 Comparison between FMT2 vs Post-FMT2

#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant3<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
  column_to_rownames("newID")

transplant3_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

full_counts<-temp_genome_counts %>%
    t()%>%
    as.data.frame()%>%
    rownames_to_column("Tube_code")%>%
    full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant3_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))

6.3.5.2 Number of samples used

[1] 49
beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)
#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]
6.3.5.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 1.180473 0.0855095 5.155229 0.035
time_point 1 0.860906 0.0623612 3.759652 0.001
type 1 1.459433 0.1057165 6.373471 0.002
Residual 45 10.304350 0.7464128 NA NA
Total 48 13.805162 1.0000000 NA NA
                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.4169018 5.739828 0.15622903   0.001      0.003   *
2   Control vs Hot_control  1 2.0940966 8.509112 0.21005427   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3004618 1.265034 0.04179854   0.141      0.423    
6.3.5.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 1.2800927 0.0939787 6.068484 0.007
time_point 1 0.9350566 0.0686477 4.432785 0.001
type 1 1.9135997 0.1404879 9.071725 0.002
Residual 45 9.4923500 0.6968858 NA NA
Total 48 13.6210990 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.8758788  8.282671 0.21084796   0.001      0.003   *
2   Control vs Hot_control  1 2.4396317 10.635546 0.24945256   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3158428  1.394345 0.04587515   0.131      0.393    
6.3.5.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1400466 0.0952654 5.873615 0.055
time_point 1 0.1138047 0.0774145 4.773017 0.001
type 1 0.1432667 0.0974558 6.008665 0.003
Residual 45 1.0729504 0.7298643 NA NA
Total 48 1.4700683 1.0000000 NA NA
                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.14387705 5.735321 0.15612552   0.001      0.003   *
2   Control vs Hot_control  1 0.22715701 9.044894 0.22036587   0.001      0.003   *
3 Treatment vs Hot_control  1 0.04648319 1.704277 0.05550617   0.134      0.402    
6.3.5.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 0.0092808 0.0077189 0.3741811 0.367
time_point 1 -0.0061674 -0.0051295 -0.2486581 0.882
type 1 0.0831052 0.0691191 3.3506286 0.257
Residual 45 1.1161295 0.9282915 NA NA
Total 48 1.2023481 1.0000000 NA NA
                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Control vs Treatment  1  0.078539743  4.59293783  0.129040706   0.066      0.198    
2   Control vs Hot_control  1  0.052468954  2.13675422  0.062593948   0.176      0.528    
3 Treatment vs Hot_control  1 -0.002340352 -0.07432315 -0.002569452   0.859      1.000    
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))
p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

6.3.5.3 Comparison between the different experimental time points (Acclimation vs Transplant samples)

The estimated time for calculating the 5151 pairwise combinations is 25 seconds.
ggarrange( p1, p2, p3, ncol=3, nrow=1, common.legend = TRUE, legend="right")

6.3.5.4 Comparison of acclimation samples to transplant samples

transplant7<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1")%>%
  column_to_rownames("newID")

transplant7_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1")

transplant7_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "1_Acclimation"| time_point == "3_Transplant1") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

transplant7_counts <- transplant7_counts[, !names(transplant7_counts) %in% c("AD45 _ LI1_2nd_2", "AD48 _ LI1_2nd_6")]

identical(sort(colnames(transplant7_counts)),sort(as.character(rownames(transplant7))))
[1] TRUE

6.3.5.5 Number of samples used

[1] 73
beta_div_richness_transplant7<-hillpair(data=transplant7_counts, q=0)
beta_div_neutral_transplant7<-hillpair(data=transplant7_counts, q=1)
beta_div_phylo_transplant7<-hillpair(data=transplant7_counts, q=1, tree=genome_tree)
beta_div_func_transplant7<-hillpair(data=transplant7_counts, q=1, dist=dist)
#Arrange of metadata dataframe
transplant7_arrange<-transplant7[labels(beta_div_neutral_transplant7$S),]
6.3.5.5.1 Richness
Df SumOfSqs R2 F Pr(>F)
Population 1 2.3108184 0.1096208 9.6707812 0.089
time_point 2 1.1082036 0.0525710 2.3189174 0.001
type 1 1.4676332 0.0696217 6.1420489 0.034
Population:time_point 2 0.4228641 0.0200599 0.8848438 0.487
Residual 66 15.7705995 0.7481267 NA NA
Total 72 21.0801189 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.3922559  4.952933 0.09533501   0.001      0.003   *
2   Control vs Hot_control  1 3.4162446 15.301454 0.23796436   0.001      0.003   *
3 Treatment vs Hot_control  1 0.6059873  2.514792 0.05406435   0.004      0.012   .
6.3.5.5.2 Neutral
Df SumOfSqs R2 F Pr(>F)
Population 1 2.641095 0.1266245 12.453113 0.038
time_point 2 1.394944 0.0668791 3.288673 0.001
type 1 1.523725 0.0730534 7.184566 0.031
time_point:type 4 1.724614 0.0826848 2.032946 0.014
Residual 64 13.573319 0.6507583 NA NA
Total 72 20.857698 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.4294917  5.324075 0.10175193   0.001      0.003   *
2   Control vs Hot_control  1 3.7896925 18.183717 0.27065661   0.001      0.003   *
3 Treatment vs Hot_control  1 0.7355091  3.012483 0.06407837   0.002      0.006   *
6.3.5.5.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
Population 1 0.3001462 0.0764583 6.402286 0.334
time_point 2 0.2984483 0.0760258 3.183035 0.003
type 1 0.1391084 0.0354361 2.967261 0.632
Residual 68 3.1879143 0.8120798 NA NA
Total 72 3.9256172 1.0000000 NA NA
                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.1229884  2.086117 0.04249913   0.089      0.267    
2   Control vs Hot_control  1 0.3936331 10.002468 0.16952626   0.001      0.003   *
3 Treatment vs Hot_control  1 0.1040675  1.985916 0.04318530   0.095      0.285    
6.3.5.5.4 Functional
Df SumOfSqs R2 F Pr(>F)
Population 1 0.0968624 0.0459940 4.124655 0.408
time_point 2 0.1660358 0.0788403 3.535121 0.040
type 1 0.2461830 0.1168973 10.483121 0.175
Residual 68 1.5968952 0.7582683 NA NA
Total 72 2.1059764 1.0000000 NA NA
                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Control vs Treatment  1 2.206030e-01 7.529285059 1.380778e-01   0.006      0.018   .
2   Control vs Hot_control  1 2.423897e-01 8.784943340 1.520282e-01   0.004      0.012   .
3 Treatment vs Hot_control  1 7.721776e-05 0.004006723 9.105359e-05   0.691      1.000    
beta_richness_nmds_transplant7 <- beta_div_richness_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant7 <- beta_div_neutral_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant7 <- beta_div_phylo_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant7 <- beta_div_func_transplant7$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(transplant7_nmds, by = join_by(sample == newID))
p0<-beta_richness_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p1<-beta_neutral_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p2<-beta_phylo_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
  theme_classic() +
  theme(legend.position="none")

p3<-beta_func_nmds_transplant7 %>%
  group_by(type) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_shape_manual(name="time_point",
                     breaks=c("1_Acclimation", "3_Transplant1", "4_Transplant2"),
                     labels=c("Acclimation", "Transplant", "Transplant"),
                     values=c("circle","square","square")) +
  geom_point(size=2) +
  geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
  labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
  theme_classic()+
  theme(legend.position="none")

6.3.5.6 Comparison between Acclimation vs Post-FMT1

post3 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post3))]
identical(sort(colnames(post3.counts)),sort(as.character(rownames(post3))))

post3_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

6.3.5.7 Number of samples used

[1] 53
beta_div_richness_post3<-hillpair(data=post3.counts, q=0)
beta_div_neutral_post3<-hillpair(data=post3.counts, q=1)
beta_div_phylo_post3<-hillpair(data=post3.counts, q=1, tree=genome_tree)
beta_div_func_post3<-hillpair(data=post3.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post3_arrange<-post3[labels(beta_div_neutral_post3$S),]
post3_arrange$type_time <- interaction(post3_arrange$type, post3_arrange$time_point)
6.3.5.7.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     2 0.099607 0.049803 9.5441    999  0.001 ***
Residuals 50 0.260911 0.005218                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Control Hot_control Treatment
Control                 0.00200000     0.886
Hot_control 0.00102653                 0.001
Treatment   0.88832670  0.00010131          
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.2059071 0.0649048 3.913237 0.001
Population 1 1.7615474 0.0948107 5.716321 0.001
time_point:Population 1 0.5122847 0.0275724 1.662393 0.010
Residual 49 15.0998916 0.8127121 NA NA
Total 52 18.5796308 1.0000000 NA NA
                                                  pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1 0.3620815 1.052109 0.06169963   0.340      1.000    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1 1.2800877 4.605444 0.22350616   0.001      0.015   .
3          Control.1_Acclimation vs Control.5_Post-FMT1  1 0.6845657 1.998114 0.11101796   0.003      0.045   .
4        Control.1_Acclimation vs Treatment.5_Post-FMT1  1 0.8437461 2.499232 0.14281954   0.001      0.015   .
5      Control.1_Acclimation vs Hot_control.5_Post-FMT1  1 1.1208022 3.568670 0.18236649   0.001      0.015   .
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1 1.3606630 5.087152 0.24124415   0.001      0.015   .
7        Treatment.1_Acclimation vs Control.5_Post-FMT1  1 0.7216200 2.172734 0.11956009   0.001      0.015   .
8      Treatment.1_Acclimation vs Treatment.5_Post-FMT1  1 0.9551308 2.926054 0.16322910   0.001      0.015   .
9    Treatment.1_Acclimation vs Hot_control.5_Post-FMT1  1 1.2263345 4.039487 0.20157637   0.001      0.015   .
10     Hot_control.1_Acclimation vs Control.5_Post-FMT1  1 1.4319792 5.384836 0.25180628   0.001      0.015   .
11   Hot_control.1_Acclimation vs Treatment.5_Post-FMT1  1 0.8172413 3.194690 0.17558364   0.001      0.015   .
12 Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1  1 0.5796135 2.441615 0.13239702   0.001      0.015   .
13         Control.5_Post-FMT1 vs Treatment.5_Post-FMT1  1 0.5615418 1.729004 0.10335366   0.016      0.240    
14       Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1 0.8438429 2.793772 0.14865413   0.002      0.030   .
15     Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1 0.3734921 1.268929 0.07799710   0.103      1.000    
6.3.5.7.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00945 0.0094472 1.1428    999  0.276
Residuals 51 0.42161 0.0082669                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.281
5_Post-FMT1          0.2901            
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.7277808 0.0986354 6.486824 0.001
Population 1 2.0558578 0.1173647 7.718565 0.001
time_point:Population 1 0.6819354 0.0389303 2.560276 0.003
Residual 49 13.0512643 0.7450696 NA NA
Total 52 17.5168383 1.0000000 NA NA
                                                  pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1 0.2316020 0.7712905 0.04598874   0.709      1.000    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1 1.4015347 5.7562378 0.26457873   0.001      0.015   .
3          Control.1_Acclimation vs Control.5_Post-FMT1  1 0.8332162 2.9081103 0.15380227   0.001      0.015   .
4        Control.1_Acclimation vs Treatment.5_Post-FMT1  1 1.1719595 4.0685514 0.21336447   0.002      0.030   .
5      Control.1_Acclimation vs Hot_control.5_Post-FMT1  1 1.4260875 5.2413171 0.24675104   0.001      0.015   .
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1 1.6347704 6.8326887 0.29925029   0.001      0.015   .
7        Treatment.1_Acclimation vs Control.5_Post-FMT1  1 0.9517634 3.3715700 0.17404733   0.001      0.015   .
8      Treatment.1_Acclimation vs Treatment.5_Post-FMT1  1 1.3127773 4.6298256 0.23585668   0.001      0.015   .
9    Treatment.1_Acclimation vs Hot_control.5_Post-FMT1  1 1.6713369 6.2395460 0.28056085   0.001      0.015   .
10     Hot_control.1_Acclimation vs Control.5_Post-FMT1  1 1.5409781 6.8338056 0.29928456   0.001      0.015   .
11   Hot_control.1_Acclimation vs Treatment.5_Post-FMT1  1 0.9133614 4.0964534 0.21451383   0.001      0.015   .
12 Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1  1 0.6954835 3.2951234 0.17077493   0.001      0.015   .
13         Control.5_Post-FMT1 vs Treatment.5_Post-FMT1  1 0.6051778 2.2508491 0.13047758   0.011      0.165    
14       Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1 1.0528902 4.1436369 0.20570451   0.001      0.015   .
15     Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1 0.4150076 1.6372683 0.09840968   0.051      0.765    
6.3.5.7.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05132 0.051320 2.6745    999  0.089 .
Residuals 51 0.97861 0.019189                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                       0.1
5_Post-FMT1         0.10812            
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.4329638 0.1243360 7.717455 0.001
Population 1 0.2375991 0.0682323 4.235135 0.002
time_point:Population 1 0.0626513 0.0179918 1.116741 0.243
Residual 49 2.7489923 0.7894398 NA NA
Total 52 3.4822065 1.0000000 NA NA
                                                  pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1 0.04186923 0.4391642 0.02671451   0.734      1.000    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1 0.15609416 2.5546889 0.13768428   0.044      0.660    
3          Control.1_Acclimation vs Control.5_Post-FMT1  1 0.19193367 2.9749922 0.15678490   0.018      0.270    
4        Control.1_Acclimation vs Treatment.5_Post-FMT1  1 0.14627288 1.7907381 0.10665035   0.155      1.000    
5      Control.1_Acclimation vs Hot_control.5_Post-FMT1  1 0.25061348 3.6146185 0.18428187   0.014      0.210    
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1 0.23108846 4.0521838 0.20208192   0.002      0.030   .
7        Treatment.1_Acclimation vs Control.5_Post-FMT1  1 0.26358465 4.3608960 0.21417997   0.002      0.030   .
8      Treatment.1_Acclimation vs Treatment.5_Post-FMT1  1 0.25319427 3.2738422 0.17915456   0.035      0.525    
9    Treatment.1_Acclimation vs Hot_control.5_Post-FMT1  1 0.39050120 5.9837393 0.27218933   0.001      0.015   .
10     Hot_control.1_Acclimation vs Control.5_Post-FMT1  1 0.14203376 5.4200212 0.25303529   0.001      0.015   .
11   Hot_control.1_Acclimation vs Treatment.5_Post-FMT1  1 0.09666753 2.3682173 0.13635351   0.023      0.345    
12 Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1  1 0.09252600 2.9824958 0.15711821   0.008      0.120    
13         Control.5_Post-FMT1 vs Treatment.5_Post-FMT1  1 0.01842535 0.4144162 0.02688498   0.785      1.000    
14       Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1 0.05987967 1.7387847 0.09802164   0.129      1.000    
15     Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1 0.03212966 0.6477782 0.04139746   0.698      1.000    
6.3.5.7.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00541 0.0054137 0.2021    999  0.659
Residuals 51 1.36615 0.0267873                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.655
5_Post-FMT1         0.65494            
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.0571984 0.0294799 1.5128525 0.174
Population 1 0.0239890 0.0123639 0.6344904 0.873
time_point:Population 1 0.0064542 0.0033265 0.1707088 0.630
Residual 49 1.8526072 0.9548297 NA NA
Total 52 1.9402488 1.0000000 NA NA
                                                  pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1  0.090450145  1.65575459  0.093779882   0.204       1.00    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1  0.085769225  1.63605364  0.092767559   0.205       1.00    
3          Control.1_Acclimation vs Control.5_Post-FMT1  1  0.030318564  0.53607587  0.032418566   0.526       1.00    
4        Control.1_Acclimation vs Treatment.5_Post-FMT1  1  0.236457683  4.06299320  0.213135113   0.038       0.57    
5      Control.1_Acclimation vs Hot_control.5_Post-FMT1  1  0.135603602  2.22854385  0.122255726   0.168       1.00    
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1  0.001373886  0.07238159  0.004503476   0.638       1.00    
7        Treatment.1_Acclimation vs Control.5_Post-FMT1  1  0.002173211  0.09402475  0.005842215   0.610       1.00    
8      Treatment.1_Acclimation vs Treatment.5_Post-FMT1  1  0.059119657  2.62461793  0.148917721   0.174       1.00    
9    Treatment.1_Acclimation vs Hot_control.5_Post-FMT1  1  0.010911935  0.39816986  0.024281360   0.500       1.00    
10     Hot_control.1_Acclimation vs Control.5_Post-FMT1  1 -0.002303582 -0.11016709 -0.006933181   0.747       1.00    
11   Hot_control.1_Acclimation vs Treatment.5_Post-FMT1  1  0.058908140  2.91987617  0.162940644   0.150       1.00    
12 Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1  1  0.005904278  0.23427876  0.014431116   0.534       1.00    
13         Control.5_Post-FMT1 vs Treatment.5_Post-FMT1  1  0.116146557  4.72479132  0.239535681   0.094       1.00    
14       Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1  0.050009298  1.70482607  0.096291602   0.209       1.00    
15     Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1  1  0.012358590  0.42381202  0.027477774   0.461       1.00    
beta_richness_nmds_post3 <- beta_div_richness_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- beta_div_neutral_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- beta_div_phylo_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- beta_div_func_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

6.3.5.8 Comparison between Acclimation vs Post-FMT2

post4 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post4.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post4))]
identical(sort(colnames(post4.counts)),sort(as.character(rownames(post4))))

post4_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

6.3.5.9 Number of samples used

[1] 54
beta_div_richness_post4<-hillpair(data=post4.counts, q=0)
beta_div_neutral_post4<-hillpair(data=post4.counts, q=1)
beta_div_phylo_post4<-hillpair(data=post4.counts, q=1, tree=genome_tree)
beta_div_func_post4<-hillpair(data=post4.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post4_arrange<-post4[labels(beta_div_neutral_post4$S),]
post4_arrange$type_time <- interaction(post4_arrange$type, post4_arrange$time_point)
6.3.5.9.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.06809 0.034047 3.8471    999   0.04 *
Residuals 51 0.45135 0.008850                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0340000     0.890
Hot_control 0.0349385                 0.005
Treatment   0.8855174   0.0047257          
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.8124061 0.0462232 2.847438 0.001
Population 1 2.0491994 0.1165926 7.182331 0.001
time_point:Population 1 0.4485668 0.0255219 1.572202 0.005
Residual 50 14.2655595 0.8116623 NA NA
Total 53 17.5757317 1.0000000 NA NA
                                                  pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1 0.3620815 1.052109 0.06169963   0.348      1.000    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1 1.2800877 4.605444 0.22350616   0.001      0.015   .
3        Control.1_Acclimation vs Treatment.6_Post-FMT2  1 0.8430295 2.845779 0.15100353   0.001      0.015   .
4          Control.1_Acclimation vs Control.6_Post-FMT2  1 0.5232174 1.683240 0.09518843   0.021      0.315    
5      Control.1_Acclimation vs Hot_control.6_Post-FMT2  1 1.1217138 3.634271 0.18509835   0.001      0.015   .
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1 1.3606630 5.087152 0.24124415   0.001      0.015   .
7      Treatment.1_Acclimation vs Treatment.6_Post-FMT2  1 0.9130048 3.195028 0.16645080   0.001      0.015   .
8        Treatment.1_Acclimation vs Control.6_Post-FMT2  1 0.5959230 1.984036 0.11032208   0.003      0.045   .
9    Treatment.1_Acclimation vs Hot_control.6_Post-FMT2  1 1.2747787 4.275366 0.21086503   0.001      0.015   .
10   Hot_control.1_Acclimation vs Treatment.6_Post-FMT2  1 0.6397330 2.913695 0.15405213   0.001      0.015   .
11     Hot_control.1_Acclimation vs Control.6_Post-FMT2  1 1.4575447 6.224524 0.28007456   0.001      0.015   .
12 Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2  1 0.3276169 1.412318 0.08111028   0.037      0.555    
13         Treatment.6_Post-FMT2 vs Control.6_Post-FMT2  1 0.6463814 2.560441 0.13795154   0.001      0.015   .
14     Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 0.4796256 1.916520 0.10696943   0.001      0.015   .
15       Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 1.1305044 4.268317 0.21059061   0.001      0.015   .
6.3.5.9.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01544 0.0154447 2.0972    999  0.151
Residuals 52 0.38294 0.0073643                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.154
6_Post-FMT2         0.15357            
Df SumOfSqs R2 F Pr(>F)
time_point 1 1.0151664 0.0602602 3.909554 0.001
Population 1 2.2827471 0.1355037 8.791191 0.001
time_point:Population 1 0.5653146 0.0335570 2.177109 0.002
Residual 50 12.9831505 0.7706790 NA NA
Total 53 16.8463787 1.0000000 NA NA
                                                  pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1 0.2316020 0.7712905 0.04598874   0.731      1.000    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1 1.4015347 5.7562378 0.26457873   0.001      0.015   .
3        Control.1_Acclimation vs Treatment.6_Post-FMT2  1 1.1746426 4.5564741 0.22165640   0.001      0.015   .
4          Control.1_Acclimation vs Control.6_Post-FMT2  1 0.5286441 1.9819408 0.11021840   0.001      0.015   .
5      Control.1_Acclimation vs Hot_control.6_Post-FMT2  1 1.3443224 4.9104417 0.23483204   0.001      0.015   .
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1 1.6347704 6.8326887 0.29925029   0.001      0.015   .
7      Treatment.1_Acclimation vs Treatment.6_Post-FMT2  1 1.3540292 5.3398081 0.25022756   0.001      0.015   .
8        Treatment.1_Acclimation vs Control.6_Post-FMT2  1 0.6311089 2.4041625 0.13063146   0.003      0.045   .
9    Treatment.1_Acclimation vs Hot_control.6_Post-FMT2  1 1.6125755 5.9825981 0.27215155   0.001      0.015   .
10   Hot_control.1_Acclimation vs Treatment.6_Post-FMT2  1 0.6202327 3.1519868 0.16457754   0.001      0.015   .
11     Hot_control.1_Acclimation vs Control.6_Post-FMT2  1 1.5701179 7.6327037 0.32297209   0.001      0.015   .
12 Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2  1 0.3634438 1.7083388 0.09647087   0.033      0.495    
13         Treatment.6_Post-FMT2 vs Control.6_Post-FMT2  1 1.0227481 4.6483346 0.22511910   0.001      0.015   .
14     Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 0.5010202 2.2065321 0.12119453   0.002      0.030   .
15       Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 1.3619424 5.7710313 0.26507845   0.001      0.015   .
6.3.5.9.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.06978 0.069777 5.0345    999  0.026 *
Residuals 52 0.72071 0.013860                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.023
6_Post-FMT2        0.029131            
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.4131659 0.1296152 8.500864 0.001
Population 1 0.2372445 0.0744265 4.881291 0.001
time_point:Population 1 0.1070826 0.0335931 2.203218 0.034
Residual 50 2.4301410 0.7623651 NA NA
Total 53 3.1876340 1.0000000 NA NA
                                                  pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1 0.04186923 0.4391642 0.02671451   0.741      1.000    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1 0.15609416 2.5546889 0.13768428   0.025      0.375    
3        Control.1_Acclimation vs Treatment.6_Post-FMT2  1 0.26322331 4.3060281 0.21205664   0.003      0.045   .
4          Control.1_Acclimation vs Control.6_Post-FMT2  1 0.16047895 2.5405742 0.13702781   0.037      0.555    
5      Control.1_Acclimation vs Hot_control.6_Post-FMT2  1 0.25529510 4.0109138 0.20043631   0.002      0.030   .
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1 0.23108846 4.0521838 0.20208192   0.001      0.015   .
7      Treatment.1_Acclimation vs Treatment.6_Post-FMT2  1 0.36496892 6.3966666 0.28560797   0.001      0.015   .
8        Treatment.1_Acclimation vs Control.6_Post-FMT2  1 0.22628210 3.8292220 0.19311005   0.016      0.240    
9    Treatment.1_Acclimation vs Hot_control.6_Post-FMT2  1 0.34830814 5.8463335 0.26761166   0.001      0.015   .
10   Hot_control.1_Acclimation vs Treatment.6_Post-FMT2  1 0.10002871 4.3836237 0.21505615   0.001      0.015   .
11     Hot_control.1_Acclimation vs Control.6_Post-FMT2  1 0.12577510 5.0601287 0.24027055   0.001      0.015   .
12 Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2  1 0.06334378 2.4997737 0.13512455   0.016      0.240    
13         Treatment.6_Post-FMT2 vs Control.6_Post-FMT2  1 0.05927454 2.3820253 0.12958449   0.030      0.450    
14     Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 0.06906280 2.7224602 0.14541146   0.002      0.030   .
15       Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 0.11081709 4.0436561 0.20174244   0.001      0.015   .
6.3.5.9.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00538 0.0053768 0.1915    999  0.675
Residuals 52 1.46001 0.0280772                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                      0.68
6_Post-FMT2         0.66348            
Df SumOfSqs R2 F Pr(>F)
time_point 1 0.0451272 0.0242533 1.2655291 0.241
Population 1 0.0004181 0.0002247 0.0117261 0.296
time_point:Population 1 0.0321822 0.0172961 0.9025046 0.304
Residual 50 1.7829366 0.9582260 NA NA
Total 53 1.8606640 1.0000000 NA NA
                                                  pairs Df     SumsOfSqs     F.Model            R2 p.value p.adjusted sig
1      Control.1_Acclimation vs Treatment.1_Acclimation  1  9.045014e-02  1.65575459  0.0937798825   0.209          1    
2    Control.1_Acclimation vs Hot_control.1_Acclimation  1  8.576922e-02  1.63605364  0.0927675587   0.224          1    
3        Control.1_Acclimation vs Treatment.6_Post-FMT2  1  1.171666e-01  2.16526354  0.1191980254   0.165          1    
4          Control.1_Acclimation vs Control.6_Post-FMT2  1  1.134013e-01  2.16495643  0.1191831338   0.191          1    
5      Control.1_Acclimation vs Hot_control.6_Post-FMT2  1  6.533744e-02  0.94814395  0.0559438221   0.264          1    
6  Treatment.1_Acclimation vs Hot_control.1_Acclimation  1  1.373886e-03  0.07238159  0.0045034760   0.634          1    
7      Treatment.1_Acclimation vs Treatment.6_Post-FMT2  1 -4.929805e-04 -0.02385162 -0.0014929516   0.716          1    
8        Treatment.1_Acclimation vs Control.6_Post-FMT2  1  2.443617e-03  0.12903840  0.0080003779   0.552          1    
9    Treatment.1_Acclimation vs Hot_control.6_Post-FMT2  1  6.968380e-03  0.19647180  0.0121305310   0.583          1    
10   Hot_control.1_Acclimation vs Treatment.6_Post-FMT2  1  3.866607e-04  0.02093980  0.0013070267   0.688          1    
11     Hot_control.1_Acclimation vs Control.6_Post-FMT2  1 -5.633463e-05 -0.00336651 -0.0002104511   0.676          1    
12 Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2  1  2.011448e-03  0.06046867  0.0037650628   0.750          1    
13         Treatment.6_Post-FMT2 vs Control.6_Post-FMT2  1 -8.527330e-03 -0.46290555 -0.0297935723   0.838          1    
14     Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1 -1.648721e-03 -0.04717131 -0.0029569243   0.905          1    
15       Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2  1  4.367477e-03  0.13147026  0.0081499244   0.680          1    
beta_richness_nmds_post4 <- beta_div_richness_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- beta_div_neutral_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- beta_div_phylo_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- beta_div_func_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

6.3.6 6. Are there differences between the control and the treatment group?

6.3.6.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

6.3.6.2 Number of samples used

[1] 26
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]
6.3.6.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.017675 0.0088373 2.3825    999  0.088 .
Residuals 23 0.085312 0.0037092                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0080000     0.625
Hot_control 0.0068795                 0.213
Treatment   0.6248469   0.2084296          
Df SumOfSqs R2 F Pr(>F)
species 1 0.6340254 0.0768024 2.065607 0.004
type 1 0.5615418 0.0680222 1.829461 0.010
Residual 23 7.0597099 0.8551754 NA NA
Total 25 8.2552771 1.0000000 NA NA
                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.022      0.066    
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.002      0.006   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.112      0.336    
6.3.6.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.566
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.22600     0.967
Hot_control 0.21166                 0.438
Treatment   0.95468     0.43604          
Df SumOfSqs R2 F Pr(>F)
species 1 0.7907904 0.1076445 3.056657 0.001
type 1 0.6051778 0.0823784 2.339205 0.006
Residual 23 5.9503501 0.8099772 NA NA
Total 25 7.3463184 1.0000000 NA NA
                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.018      0.054    
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.001      0.003   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.054      0.162    
6.3.6.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.899
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92100     0.665
Hot_control 0.91505                 0.758
Treatment   0.63312     0.73046          
Df SumOfSqs R2 F Pr(>F)
species 1 0.0560850 0.0531376 1.3149967 0.269
type 1 0.0184254 0.0174571 0.4320099 0.814
Residual 23 0.9809570 0.9294053 NA NA
Total 25 1.0554673 1.0000000 NA NA
                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.762      1.000    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.106      0.318    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.694      1.000    
6.3.6.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.00400 0.0020014 0.145    999  0.859
Residuals 23 0.31753 0.0138057                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.58000     0.730
Hot_control 0.59817                 0.836
Treatment   0.75141     0.83718          
Df SumOfSqs R2 F Pr(>F)
species 1 0.0024979 0.0033024 0.0900845 0.629
type 1 0.1161466 0.1535542 4.1887855 0.058
Residual 23 0.6377435 0.8431434 NA NA
Total 25 0.7563879 1.0000000 NA NA
                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11614656 4.724791 0.23953568   0.084      0.252    
2   Control vs Hot_control  1 0.05000930 1.704826 0.09629160   0.224      0.672    
3 Treatment vs Hot_control  1 0.01235859 0.423812 0.02747777   0.495      1.000    
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.6.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.6.4 Number of samples used

[1] 27
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]
6.3.6.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.841
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.69800     0.824
Hot_control 0.67789                 0.638
Treatment   0.79246     0.59820          
Df SumOfSqs R2 F Pr(>F)
type 2 1.504341 0.1967776 2.939822 0.001
Residual 24 6.140538 0.8032224 NA NA
Total 26 7.644879 1.0000000 NA NA
                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 0.6463814 2.560441 0.1379515   0.001      0.003   *
2 Treatment vs Hot_control  1 0.4796256 1.916520 0.1069694   0.001      0.003   *
3   Control vs Hot_control  1 1.1305044 4.268317 0.2105906   0.001      0.003   *
6.3.6.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999   0.46
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.45700     0.671
Hot_control 0.44675                 0.238
Treatment   0.65989     0.25095          
Df SumOfSqs R2 F Pr(>F)
type 2 1.923807 0.2603795 4.224537 0.001
Residual 24 5.464666 0.7396205 NA NA
Total 26 7.388473 1.0000000 NA NA
                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 1.0227481 4.648335 0.2251191   0.001      0.003   *
2 Treatment vs Hot_control  1 0.5010202 2.206532 0.1211945   0.002      0.006   *
3   Control vs Hot_control  1 1.3619424 5.771031 0.2650785   0.001      0.003   *
6.3.6.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.947
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93400     0.839
Hot_control 0.93765                 0.774
Treatment   0.83933     0.76015          
Df SumOfSqs R2 F Pr(>F)
type 2 0.1594363 0.2042241 3.079623 0.001
Residual 24 0.6212564 0.7957759 NA NA
Total 26 0.7806927 1.0000000 NA NA
                     pairs Df  SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 0.05927454 2.382025 0.1295845   0.025      0.075    
2 Treatment vs Hot_control  1 0.06906280 2.722460 0.1454115   0.010      0.030   .
3   Control vs Hot_control  1 0.11081709 4.043656 0.2017424   0.005      0.015   .
6.3.6.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01126 0.0056302 0.2861    999  0.796
Residuals 24 0.47233 0.0196806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.51200     0.660
Hot_control 0.48255                 0.816
Treatment   0.60116     0.75643          
Df SumOfSqs R2 F Pr(>F)
type 2 -0.0038724 -0.0056213 -0.0670788 0.915
Residual 24 0.6927468 1.0056213 NA NA
Total 26 0.6888744 1.0000000 NA NA
                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Treatment vs Control  1 -0.008527330 -0.46290555 -0.029793572   0.877          1    
2 Treatment vs Hot_control  1 -0.001648721 -0.04717131 -0.002956924   0.919          1    
3   Control vs Hot_control  1  0.004367477  0.13147026  0.008149924   0.723          1    
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")